Ab-initio path intégral techniques for molécules 
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Path intégral Monte Carlo with Green's function analysis allows the sampling of quantum mechan- 
ical properties of molécules at finite température. While a high-precision computation of the energy 
of the Born-Oppenheimer surface from path intégral Monte Carlo is quite costly, we can extract 
many properties without explicitly calculating the electronic énergies. We demonstrate how physi- 
cally relevant quantifies, such as bond-length, vibrational spectra, and polarizabilities of molécules 
may be sampled directly from the path intégral simulation using Matsubura (température) Green's 
functions (imaginary-time corrélation functions). Thèse calculations on the hydrogen molécule (H2) 
are a proof-of-concept, designed to motivate new work on fixed-node path- intégral calculations for 
molécules. 



I. INTRODUCTION 

The quantum Monte Carlo (QMC) method has been 
extensively used as a powcrful technique for calculating 
some physica]_properties of many-body Systems ŒL such 
as molécules 0, Q , solids [3] , and nanostructures . For 
molécules, it has been reported that the diffusion QMC 
method can achieve near chemical accuracy within the 
fixed-node approximation 0] . 

While variational and diffusion QMC algorithms usu- 
ally sample only électrons at zéro température, path in- 
tégral Monte Carlo (PIMC) can sample both the ions 
and électrons simultaneously at finite température. This 
is a promising approach to calculate thermal effects in 
the molecular Systems. Additionally, PIMC includes the 
zero-point motion of the ions, yeilding the correct quan- 
tum statistical treatment of the motion of light ions at 
ail températures. 

We note that another approach, coupled electron-ionic 
Monte Carlo (CEIMC) [1, Q, has been demonstrated in 
which électrons and ions are sampled in two independent 
but coupled Monte Carlo algorithms. For the molécule 
we study in this paper, a comparable CEIMC simulation 
would have a ground-state variational or diffusion QMC 
simulation of the électrons that détermines the Born- 
Oppenheimer energy surface. Thèse QMC électron éner- 
gies would then be used in a path intégral simulation for 
the protons. Because the QMC estimâtes of ground-state 
electronic énergies have statistical error bars, a penalty- 
method 3 would be used to partially compensate for 
non-linear biases when thèse énergies are fed into the 
ionic path intégral Monte Carlo. Ultimately, the CEIMC 
method requires a fairly accurate détermination of total 
énergies at many points on the Born-Oppenheimer sur- 
face for many différent configurations of the ion positions. 

In this paper we addrcss the foUowing issue: by treat- 
ing électrons and ions in the same path intégral sim- 
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ulation (avoiding the explicit calculation of the Born- 
Oppenheimer energy surface), can we extract properties 
of a molécule, such as bond-lengths, vibrational frequen- 
cies, and polarizability? Thèse properties are often cal- 
culated directly from the Born-Openheimer surface, i.e. 
bond-lengths and vibrational frequencies come from lo- 
cal minima and curvature of the energy surface. Here 
we treat thèse properties as expectation values of the 
electron-ion thermal density matrix. That is, we do not 
make the Born-Oppenheimer approximation, instead, we 
sample the électron and ion coordinates in a single path 
intégral. The calculations presented here are a proof- 
of-concept, demonstrating how to extract physical prop- 
erties from the PIMC simulations. The question of the 
relative efficiencies of this approach and CEIMC will have 
to be left to future calculations on larger Systems, most 
likely requiring some kind of a path-integral fixed-node 
approximation [10]. 

In PIMC simulations we can coUect varions corrélation 
function of the System in imaginary time. Thèse func- 
tions, known as finite température Matsubura Green's 
functions, can be transformed from the imaginary-time 
domain into the imaginary-frequencics through the fast 
Fourier transform (FFT). Since we want to calculate 
the physical properties of the Systems, the corrélation 
function at real frequency need to be obtained from the 
imaginary-frequency function. It is possible through the 
analytic continuation techniques to obtain real-time cor- 
relation function from imaginary-time corrélation func- 
tion, since a time-ordered corrélation function is an ana- 
lytic function of the time variable in the complex plane 
[11.] . In othcr words, the time-correlation function calcu- 
lated along the imaginary axis can be uniquely analyti- 
cally continued to the real-time axis. 

For the case of the simple and closed corrélation func- 
tion, the analytic continuation of imaginary-frequency 
corrélation function is simply the replacement of icun with 
Lu + iri. But for realistic corrélation functions in molécules 
it is hard to find the closed form of imaginary-frequency 
corrélation function. Therefore analytic continuations 
need to be performed with numerical methods. 

However, the analytic continuation of numerically eval- 
uated corrélation functions often leads to grossly ampli- 
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fied statistical errors in the real time corrélation fonction. 
Some numerical methods have been proposed for reduc- 
ing the statistical errors, such as using Padé approxima- 
tions [HI or least-squares fitting [13], minimizinga suit- 
ably defined potential [14i] . maximum-entropy ^ISl . 
and singular value décomposition (SVD) methods 
It has been reported that the rotational constant of a 
molécule can be obtained by fitting iniaginary-time cor- 
relation function from PIMC data to analytic models of 
corrélation function [ï^ . 

In this paper we show that some properties of 
molécules can be sampled directly from path intégral sim- 
ulations using Matsubura Green's functions. For an accu- 
rate treatment of the Coulomb interactions between par- 
ticles we use a new method to tabulate of the Coulomb 
density matrix that we have developed. Using corréla- 
tion functions sampled with the accurate Coulomb prop- 
agator, we demonstrate PIMC calculations of the bond 
length, vibrational frequency, and polarizability of a hy- 
drogen molécule. This calculations are important bench- 
marks for future calculations on more complicated Sys- 
tems. In particular, we have chosen a simple System with- 
out a fermion sign-problem so that we can separate the 
analysis of corrélation functions from questions about the 
fixed-node approximation. 



II. MONTE CARLO SAMPLING OF PATH 
INTEGRALS 

A. Monte Carlo Sampling 

In thermal equilibrium, the average value of any phys- 
ical quantity can be calculated from the density matrix. 
In the coordinate représentation, the thermal density ma- 
trix of N particles is defined by 



p(R,R';/3) = i(R|e-^'«|R'), 



(1) 



where /3 = l/ksT is the inverse température, R = 
(ri,--- ,rAr) are the particle coordinates, is the po- 
sition of the i-th particle, and Z — J p(R, R; /3)(iR is the 
partition function. 

The average value of any physical quantity O can be 
written 



(O) 



dRdR'/9(R,R';/3)(R|C|R')- 



(2) 



In the path-integral formula, in N slices, we expand the 
density matrix [19| , 



p(Ro,R7v;/?) = J dRidR2---dRjv-i^ 

(Rn-l — Rn)^ 



-3/2 



X exp 



N 

E 

ra=l 



ArF(R„) 



where m is the mass of the particle and the time step 
is At — f3/N. In our molécule simulations Ar = 
0.01 Ha"^. We estimate this high dimensional intégral 
with Metropolis Monte Carlo algorithm, with an efficient 
multilcvel sampling method. For a detailed discussion of 
PIMC methods, see référence [201 • 



B. Accurate Coulomb Action 

The primitive approximation to the path intégral, 
Eq. ([3]), is not correct for attractive Coulomb interac- 
tions. Instead, we use an improved quantum action 
U (R, R') in place of the bare Coulomb potential V^(R). 
The use of improved actions is discussed in référence |2Çi] . 

In order to more accurately describe the Coulomb in- 
teractions we have developed a new technique for fabu- 
lation of imaginary-time Coulomb Green's function [2Ï[, 
from which we extract the Coulomb action. The Green's 
function, G(r, r') for two particles in A^-dimensional 
space satisfies the équation. 



2fi 



GW(r,r';r) 



dr 



GW(r,r';r). (4) 



The initial condition is G(^)(r,r';0) = 6^{r - r') and 
r and r' refer to the relative séparation of two particles 
with reduced mass and z is the product of the charges. 
There are several ways to numerically evaluate Green's 
function, see référence [2Ï[. But each methods has its 
particular limitation and numerical errors are often more 
than one percent. 

Since the Coulomb potential, V{r) ce i, has the 
symmetry properties, without a partial wave expansion 
we are able to simple, fast, and highly accurate nu- 
merical évaluation of G(r, r';r) using one-dimensional 
FFT's. The technique relies on Hostler's recursion re- 
lation 0, H^l , which relates the Ai"-dimensional Green's 
function to the 1-dimensional function. The use of an ac- 
curate Coulomb action is especially important in chem- 
istry. An analogy can be drawn to the use of accurate 
intégration schemes in molecular dynamics. In the cur- 
rent work, we find that quantifies, such as the molecular 
binding energy and electrical polarizability are especially 
sensitive to the quality of the tabulated Coulomb action. 



III. ESTIMATORS 

A. Energy and the Virial Estimator 

There are several estiniators to obtain the total ther- 
mal energy of a System. The most widely used energy 
estimator is the thermodynamic energy estimator, which 
is obtained by differentiating the partition function with 
respect to /3, 



(3) 



E =-l^ = 
^ Zd(3 



dU^ 
dr 



(5) 
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where Z is the partition function and Ui is the action. 
In other words, the thermal energy is an average value of 
the imaginary-time derivative of the action. 

For the calculation of molécules we used the virial en- 
ergy estimator. The virial estimator of the energy is 



Ev = 



dr 



1 



(6) 



where is the déviation of the î-th particle from its 
average position, and Fi is a generalization of the classical 
force, 



F, 



1. 



(7) 



For more detailed dérivation, see références [20| and [24 



B. Static polarization 



The gênerai description of the static polarization of 
molécules is P = a^ti^E, where E is the applied elec- 
tric field, a^^, is the static polarizability tensor of the 
molécule, and /i and v dénote vector directions. For di- 
atomic molécules, symmetry reduces the polarizability 
tensor to reduced into two components, the parallel (a||) 
and perpendicular polarizability. 

Using the polarization estimator 



a 



fil- 



(8) 



we calculate the static polarizability of the molécules 
given the electric field. Here X^{i) is the position of 
the i-th particle in the direction of fi, and Ei, is the mag- 
nitude of the the electric field in the v direction. 



where u;„ = 2n7r//3 and n = 0, ±1, ±2, ■ • ■ . In the linear 
response theory, the linear response of the System to a 
small perturbation can be written in terms of the Green's 
function. 

As mentioned in the introduction section, through ana- 
lytical continuation, Matsubura Green's function of com- 
plex frequencies, itOn, can be analytically continued to 
real, time-ordered corrélation functions, which détermine 
physical properties of the dynamic System. However, 
in numerical analytic continuation even small noise on 
imaginary-time corrélation can make big statistical er- 
rors in real time dynamics of the System |ï^. Instead 
of that, we analyze directly the imaginary-time corréla- 
tion function calculated by path intégral Monte Carlo 
method in order to obtain some properties of molécules. 
To calculate bond-length and vibrational frequency, we 
use displacement-displacement imaginary-time corréla- 
tion functions and their Fourier transforms. For polariz- 
ability of molécules we use polarization-polarization cor- 
relation function. 



1. Bond-length and vibrational frequency 

For calculating of the bond-length and vibrational 
frequencies of molécules we use the displacement- 
displacement corrélation function. The displacement cor- 
relation function of a harmonie oscillator, with V{x) — 
imc^sHo^c^, in the imaginary-time domain at the inverse 
température P can be written 



G(t) = -(r,x(r)x(0));3, 



(11) 



where x{t) refers to the displacement of the oscillator 
from equilibrium. For r > 0, this corrélation function is 

m, 



C. The Matsubura Green's function method 

The Matsubura Green's function method is a very use- 
ful technique for calculating the physical properties of 
many-body Systems at finite température [25| . The gên- 
erai définition of température Green's function is 



G(r,r')^-(r.A(r)B(r'))^, 



(9) 



where Tr is the time-ordering operator, r represent the 
imaginary time, and A and B are operators. The bracket 
means the thermal average value at the inverse tempér- 
ature (3 = l/ksT. If the Hamiltonian of the System is 
independent of time, Green's function becomes only a 
function of the time différence r — r', so it also can be 
written G(r) = -(r^A(T)B(O))0. 

A Fourier transform gives the frequency dependence of 
Green's function for bosonic quantifies, such as molecular 
vibrations and polarizabilites. 



P Jo 



(10) 



G(r) 



1 



2rncJsH0 



(iV + l)e- 



(12) 



where N = l/{. 



1) is the distribution function at 



the température, ksT = 1//3 and we set h= 1. 

For the hydrogen molécule, the displacemcnt x{t) rep- 
resents the déviation from the equilibrium séparation dis- 
tance between two protons. The mass, m is the reduced 
mass of two protons. At low température, we assume 
that the molécule is in harmonie motion, so we compare 
our results with the simple harmonie oscillator. 

We are able to extract the bond length from imaginary- 
time displacement-displacement corrélation function, 
Eq. (ini)- The function G(t) = -{x{t)x{Q)) can mea- 
sure the fluctuation in the quantity x. At r = 0, the 
initial value of G(0) is the equilibrium average — (a;^). 
As time goes on, the corrélation function shows the time 
decay which measure how x{t) and a;(0) are correlated 
each other. In the long time limit, the corrélation func- 
tion becomes totally uncorrelated so that the corrélation 
function can be written as, 



G(t 



-{x{t)){x{Q)) 



(13) 
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At r = 13/2, G(/3/2) is approximately {D)'^ , where D 
is the bond length of the molécule at low température. 
Since the bond length can be directly computed as a time- 
independent average, this limit is an important sanity 
check on our imaginary-time data. 

To obtain the vibrational frcqucncy of molécules we 
use imaginary-frequency corrélation function. Using 
Eq. pH)) . we can obtain the corrélation function of the 
simple harmonie oscillator in the imaginary-frequency 
space, 



harmonie oscillator is 



l/m 



1 

'1 



, ,2 



'Xjm 



(14) 



Once we have the corrélation function, we can calculate 
the vibrational frequency Wshd at iujn — 0, 



l/(/3m) 
G(zc^„ = 0) 



(15) 



at each températures (3. 

Another way to find the frequency is to directly use 
the corrélation function. Since the corrélation function 
Giiiûn) forms a Lorentzian in the imaginary axis of com- 
plex Green's function, the fuU width at half maximum 
(FWHM) of the corrélation function is equal to the fre- 
quency of the harmonie oscillator. Therefore we can com- 
pare to thèse frequencies using two methods. 



2. Polarization 

We use the polarization-polarization corrélation func- 
tion to calculate the polarizability of molécules. The 
imaginary-time corrélation function can be written 



(16) 



where P^t) — ex^r). Once we have obtained the 
imaginary-time corrélation function, by Fourier trans- 
form Eq. (IlOp we can calculated the dynamic polariz- 
ability. 

Similar to the bond-length calculation, first we calcu- 
late the polarizability for a simple harmonie oscillator. 
(Note that this oscillator is a model for our polariza- 
tion analysis and is not related to molecular vibrations.) 
From lincar responses theory the static susceptibility a, 
due to the small electric field perturbation —exE, can be 
written with the first order perturbation theory. 



a 



2e2(0|x2|0) 



(17) 



It is clear that the polarizability dépends on the confine- 
ment potential of the harmonie oscillator, o^shd- 

We can also obtain the polarizability from the 
imaginary-frequency polarization-polarization corréla- 
tion function. The corrélation function for the simple 



a(ia;„) = / dre''-"^(T,P,,(r)P,(0)) 



2/, 



(18) 



(lW„)2 - uji 



Similar to bond length estimator, we assume that the 
molécules are in the harmonie motion. Once we obtain 
a{iLUn), the polarizability a{0) is the value at iujn — 0. 
Note that for molécules the polarizations P are to be 
added for ail particles in the molécule. 



IV. RESULTS 



In 



this section we report the results of our tests 
on an H2 molécule. Simulations were performed in 
sériai and in parallel jobs with up to eight proces- 
sors using our open-source pi code available online at 
http://phy.asu.edu/shumway/codes/pi.html, We use a 
time step of Ar — 0.01 Ha^^, which leads to a discretiza- 
tion of the path intégral into 100,000 slices at the lowest 
température (when /3 = 10^ Ha-^ for T « 300 K). 



A. Born-Oppenheimer energy surface 

To check the accuracy of our discretized path intégral, 
we calculate the potential energy surface of the ground 
State H2 molécule with the internuclear séparation from 
0.5 ao to 4.5 oq. For this calculation we fixed the proton's 
position of the H2 molécule, then calculate the thermal 
energy using the virial energy estimator. Figure [T] shows 
a very good agreement with the accurate potential energy 
curve(solid line) [27]. Note that we do not use this energy 
surface for any of the subséquent calculations. We only 
présent the calculation here to verify that we can get 
the correct detailed energy surface if we so désire. Ail 
foUowing results are directly sampled from the electron- 
ion density matrix. 



B. Vibrational frequency of a H2 molécule 

By sampling the displacement-displacement corréla- 
tion function in imaginary time, we have calculated 
the vibration frequencies of a hydrogen molécule (H2) 
at three différent inverse températures, (3 = 200, 500, 
and 1000 Ha"!, corresponding to T w 1500K, 600K, 
and 300K. Figure [2] shows the corrélation function in 
imaginary-frequency domain for a hydrogen molécule. 
The solid lines are the displacement-displacement cor- 
relation functions of imaginary-frequency for the simple 
harmonie oscillator at each températures and the points 
represent our data points for the Hydrogen molécule. As 
described in Sec. IIII C 11 we extract the frequency using 
two différent approaches: (i) a fit to a simple harmonie 
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FIG. 1: Born-Oppenheimer surface of H2 molécule. Solid 
line is the accurate potential energy curve from the référence 
and data points are PIMC results using the virial energy 
estimator. 



TABLE I: Vibrational frequencies calculated from PIMC sim- 
ulations, using a fit to a simple harmonie oscillator (SHO fit) 
and a measurement of the fuU-width at half-maximum in the 
imaginary frequency response. 



/3(Ha-i) 


T(K) 


huj = El - 
SHO fit 


-Eo 

FWHM 


200 


1500 


0.01866(5) 


0.01867 


500 


600 


0.01808(5) 


0.01810 


1000 


300 


0.01786((5) 


0.01788 



oscillator and (ii) the fuU-width at half-maximum of the 
frequency response. The results are summarized in Ta- 
ble E 

The exact ground state vibrational frequency of H2 
is 0.02005 in atomic units (=4401.21 cm^i) 28]. Due 
to the anharmonic properties of the hydrogen molecule's 
potential surface, the différence between the ground state 
and the first excite state level is 0.01895 hartree (4161.14 
cm~^) [2^. The time corrélation functions measure 
this physical energy différence, not the unphysical exact 
ground state frequency. 

It is shown that the frequency is not sensitive to the 
températures. Because the thermal energy is so small, 
it just affects the rotational levels. Since we doesn't fix 
the protons in the space, there are rotational effects to 
the calculation of frequency. The analytical results (solid 
lines in Fig. [2]) are based on the harmonie oscillator with- 
out rotational motion. 



FIG. 2: The displacement-displacement corrélation functions 
of imaginary- frequency for a hydrogen molécule at three dif- 
férent températures. The lines are the corrélation functions 
of imaginary-frequency of harmonie oscillator at each tem- 
pératures and points represent our results before fitting the 
vibrational frequencyinto the corrélation function. 



the displacement-displacement corrélation functions of 
imaginary time at three températures. At room tempra- 
ture, Fig. [3] (c), we could oberve the flat in the middle of 
the imaginary time, which shows the corrleation function 
is fuUy uncorrelated. 

The exact bond length of a ground state hydrogen 
molécule is 1.401 in atomic units [2^|. With bond-length 
estimator we directly calculate two average values of the 
internuclear distance, (D) and (D^). Table [TTl shows the 
summary of two bond length averages from the estima- 
tor and the corrélation function of a hydrogen molécule 
in atomic units. It shows the two différent methods agrée 
very well. 

The thermal effects on the bond-length of a molécule 
can be calculated from the Boltzmann factor (2J -|- 
l)e~^^-' with 2 J + 1-fold degeneracy due to the rota- 
tional motion, where J is the total angular momentum. 
The energy is Ej = BJ{J+ 1), where B is the rotational 
constants of a hydrogen molécule, B = 60.853 cm~^. At 
température, fcsT = 1//3, the number of molécules Nj 
in the rotational level J can be calculated by 



Nj 



(2J+l)e-^^'^('^+i) 



(19) 



C. Bond-length of a H2 molécule 

Next, we calculate the bond length of the H2 molécule 
from imaginary-time corrélation function. Figure[3]shows 



With the number of distribution with respect to the tem- 
pérature, we can calculate the effective potential due to 
the rotational motion in explicit method using Morse po- 
tential. We observed that the bond-length has changed 
the factor of 10~^. 
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FIG. 3: The displacement-displacement corrélation function 
of imaginary time for a hydrogen molécule at three différent 
températures. The solid lines are analytical corrélation func- 
tion of simple harmonie oscillator and the points are PIMC 
results at each températures. 



TABLE II: The summary of two bond length averages from 
the bond-length estimator and the displacement-displacement 
corrélation function of a hydrogen molécule (IÎ2) by path in- 
tégral simulation (PIMC CF) in atomic units. 







Bond-lenj 


;th estimator 


PIMC CF 


/3(Ha-i) 


T(K) 


(Df 


{D^} 


at T = /J/2 at r = 


200 


1500 


2.070(4) 


2.101(5) 


2.076(1) 2.097(1) 


500 


600 


2.029(4) 


2.056(5) 


2.031(1) 2.057(1) 


1000 


300 


2.024(4) 


2.051(5) 


2.029(1) 2.055(1) 



D. Polarization 

Many-body perturbation methods are usually used for 
calculating the polarizability of molécules. In the case of 
the hydrogen molécule, the wavefunctions are very well 
known so that the variational perturbation methods of 
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FIG. 4: The polarization-polarization corrélation function 
for simple harmonie oscillators with three différent harmonie 
strength, usw = 0.5, 1, 2. At itOn = 0, the arrows show the 
polarizability for each oscillators. 



the polarizability of the molécule agrée very well with 
the expérimental results [s^, [3Ï| . With the same method 
the polarizability have been calculated of the depedence 
on rotation and vibration states [32] and of the function 
of internuclear séparation (ssj . 

We have calculated the polarizability of a hydrogen 
molécule in PIMC using two différent approaches. The 
first method is to use the polarization estimator as we 
explained in the estimator section, see Eq ([5]). By ap- 
plying a small electric field through the z-axis to a fixed 
molécule, we have calculated the polarzability at room 
température, for the parallel direction a|| = 6.12(50), and 
for perpendicular a±= 4.92(26) in atomic units. The ex- 
act value is a|| =6.380 and a^=4.577 [31]. 

Anothcr way to calculate the polarizability is to use 
the polarization corrélation function. We calculate the 
polarizability of the simple harmonie oscillator for three 
différent oscillator strength wsho = 0.5,1,2 for the test 
of the corrélation function method. As we mention in 
the estimator section, the static polarizability can be ob- 
tained at iujn = of a(ia;„). As the harmonie strength 
increase, the polarizability becomes smaller, see the équa- 
tion ^TE\i . Figure 2] shows that the values of polarizabil- 
ity for the différent harmonie strengths exactly foUow 
Eq. (US]). The polarizabilities for each oscillators are for 
t^sHO = 1, Q! = 1, for WsHO = 0.5, a — A and for luseo — 2, 
a = 0.25 respectively at icun — 0. 

We can also estimate the possible dipole transition en- 
ergy levels (husuo) h'om the values of the fuU width at half 
maximum (FWHM) of corrélation function. With the an- 
alytic continuation we can obtain retarded Green's func- 
tion (G^(w)) of real frequency (w) from the imaginary- 
frequency corrélation function (G'(iw„)). The retarded 
Green's function is also related the spectral density func- 



7 



a -2 




FIG. 5: The polarization-polarization corrélation function of a 
hydrogen atom. At iu>n = 0, the arrow shows the polarzabihty 
of a hydrogen atom, 4.51(3). The exact value is 4.5 in atomic 
units. 



tion {A{uj)) for the harmonie oscillator at > 0, 



duj' - 



1 

1/m 
{uj + iriY - 



IT] 



(20) 



, .2 ' 



At the pôles lu ~ cj^sho, the speetral function is a delta 
function and imaginary-frequency corrélation function 
form a Lorentian with FWHM of o^sho- Since the har- 
monie oscillator has one possible transition, the values of 
FWHM are exactly the same as the harmonie strengths. 
But for a atom and molécule the FWHM doesn't repre- 
sent the possible transition levels, because the levels are 
not équivalent. 

Next we test the corrélation function mcthod for cal- 
culating the polarizability of a hydrogen atom. We 
found that the polarization-polarzation corrélation func- 
tion has the same form whether the proton of the atom 
is fixed or not. Figure [5] shows the corrélation function 
at /3 = 1000 Ha~^. We have obtained the polarizability 
is 4.51(3) in atomic units, while the exact value is 4.5. 

For the hydrogen molécule, we fixed the positions of 
two protons at the z-axis, so that the rotation effect 
doesn't need to be considered for the polarization of the 
molécule and thus the corrélation function is independent 
of the température. 

Figure [S] shows the imginary-frequency corrélation 
function of a hydrogen molécule when the protons are 
fixed (a) and unfixed (b). For the perpendicular (a^) 
and parallel {ol\\) polarizability we need to calculate x,y 
and z components of the displacement of particles(two 
électrons and two protons) respectively. We have ob- 
tained the perpendicular aj^ = 4.56(17) and parallel 
«Il = 6.40(15) polarizability. 

The value of FWHM is 0.596 for the perpendicular 
(ax) and 0.534 for the parallel {ct\\) of the corrélation 
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FIG. 6: The polarizability of a hydrogen molécule at /3 = 
1000. The protons are fixed on the z-axis for the upper panel 
and unfixed at the lower panel. The x, y components and the z 
component of polarizability represent a± and ay respectively. 



function. The polarizability of a molécule is a physical 
quantity how the molécule induced by an electric field. 
For a diatomic mole possible dipole transition energy lev- 
els for {a±) is bigger than («n) because of the molécule 
structure. 



We also have calculated the corrélation function with 
unfixed protons, allowing rotational and vibrational mo- 
tion of the molécule. Since the trace of the polar- 
izability tensor is invariant, the average polarizability, 

a — — (2q;j^ + <^||)ï should be identical with différent con- 

figuration, see Fig. [5] (b) . The exact average polarzabihty 
is 5.4139 in atomic units f3Ï|. This is the same value of 
nonfixed protons. Table IHII shows the summary of the 
polarizability for the simple harmonie oscillators, the hy- 
drogen atom and the hydrogen molécule. 
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TABLE III: The summary of the polarizability of simple har- 
monie oseillators (SHO), a hydrogen atom (H), and a hydro- 
gen moleeule (H2) for the fixed and unfixed protons using path 
intégral Monte Carlo method with polarization-polarization 
corrélation function methods (PIMC-CF) in atomic units. 
The exact data of H2 are from référence l3ll [. 







PIMC-CF 


Exact 


SHO 


0.5 


4.00(6) 


4 




1 


1.00(6) 


1 




2 


0.25(6) 


0.25 


H atom 


a 


4.51(3) 


4.5 


H2 (fixed) 


ai_ 


4.56(17) 


4.577 




Q|| 


6.40(15) 


6.380 


H2 (unfixed) 


a 


5.38(18) 


5.413 



V. CONCLUSIONS 

We have demonstrated analysis techniques for path in- 
tégral Monte Carlo simulations of molécules. The aim 
of thèse calculations is to show how important physical 
properties of molécules — énergies, bond lengths, vibra- 
tional frequencies, and polarizability — can be computed 
directly from a path intégral simulation without an cx- 
plicit calculation of énergies on the Born-Opcnhcimer 
surface. This paper has focused on the analysis of Mat- 
subura Green's functions for bond length and electrical 
dipole operators. We have limited this analysis to the H2 
molécule to avoid complications from fixed-node approx- 
imations. Future calculations will need to properly treat 
fermions, most likely with a fixed-node approximation, if 
thèse techniques are to have use in practical molecular 
simulations. The results presented here are benchmarks 



that validate continued development of the PIMC tech- 
nique. 

We have two main findings: (1) displacement- 
displacement Green's functions can give accurate esti- 
mâtes of the energy splitting huj between the ground 
and first-excited vibrational states, and (2) dipole-dipole 
Green's functions provide accurate estimâtes of the elec- 
trical polarizability of atoms and molécules. Implicit in 
thèse results is the démonstration of very high accuracy 
of the discretized path intégral through the use of a care- 
fuUy tabulated coulomb action. 

It is important to the computational strategy we are 
advocating in thèse calculations: rather than using costly 
many-body techniques to calculate individual points on 
the Born-Oppenheimer surface, we are including the 
quantum and thermal average of ionic positions in the 
many-body electronic calculation. This approach may 
lead to better algorithmic efficiency for Systems that have 
significant thermal or zero-point ionic motion, includ- 
ing Systems with strong polaronic effects. The merits 
of this approach relative to other QMC techniques, such 
as CEIMC, cannot be measured until we treat larger 
molecular Systems. More complicated molécules will also 
lead to more complicated Green's functions that cncode 
information on the normal modes. Nevertheless, the 
benchmark calculations presented here are a necessary 
first step towards future calculations on more interesting 
molécules. 



Acknowledgments 

Work supported by NSF Grant No. DMR 0239819 and 
made use of facilities provided by the Ira A. Fulton High 
Performance Computing Initiative. 



[1] R. Guardiola, in Microscopic Quantum Many-Body Thé- 
ories and Their Applications, edited by J. Navarro and 
A. Poils (Springer, 1998), pp. 269-336. 

[2] P. J. Reynolds, D. M. Ceperley, B. J. Aider, and 
J. William A. Lester, J. Chem. Phys. 77, 5593 (1982). 

[3] B. L. Hammond, J. W. A. Lester, and P. J. Reynolds, 
Monte Carlo methods in Ah Initia quantum chemistry 
(World Scientific, 1994). 

[4] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Ra- 
jagopal, Rev. Mod. Phys 73, 33 (2001). 

[5] J. Shumway and D. M. Ceperley, in Handbook of Theoret- 
ical and Computational Nanotechnology, edited by M. Ri- 
eth and W. Schommers (American Scientific Publishers, 
Germany, 2006). 

[6] J. C. Grossman, J. Chem. Phys. 117, 1434 (2002). 

[7] N. A. Benedek, I. K. Snook, M. D. Towler, and R. J. 
Needs, J. Chem. Phys. 125, 104302 (2006). 

[8] D. M. Ceperley, M. Dewing, and C. Pierleoni, in Tap- 
ies in Candensed Matter Physics, edited by P. Nielaba, 



M. Mareschal, and G. Ciccotti (Springer- Verlag, 2002). 
[9] C. Pierleoni and D. M. Ceperley, ChemPhysChem 6, 
1872 (2006). 

[10] J. Shumway, in Computer Simulations Studies in Can- 
densed Matter Physics, edited by D. P. Landau, S. P. 
Lewis, and H.-B. Schiittler (Springer, Berlin, 2005), vol. 
XVII. 

[11] G. Baym and N. D. Mermin, J. Math. Phys. 2, 232 
(1961). 

[12] D. Thirumalai and B. J. Berne, J. Chem. Phys. 79, 5029 
(1983). 

[13] H. B. Schuttler and D. J. Scalapino, Phys. Rev. Lett. 55, 
1204 (1985). 

[14] M. Jarrell and G. Biham, Phys. Rev. Lett. 63, 2504 
(1989). 

[15] E. Gallicchio and B. J. Berne, J. Chem. Phys. 105, 7064 
(1996). 

[16] G. Krilov, E. Sim, and B. J.Berne, J. Chem. Phys. 114, 
1075 (2001). 



9 



[17] E. Gallicchio, S. A. Egorov, and B. J. Berne, J. Chem. 

Phys. 109, 7745 (1998). 
[18] N. Blinov, X. Song, and P.-N. Roy, J. Chem. Phys. 120, 

5916 (2004). 

[19] P. R. Feynman, Statistical Mechanics (Westview Press, 
1972). 

[20] D. M. Ccporlcy, Rev. Mod. Phys 67, 279 (1995). 

[21] D. Shin, ,]. Shumway, and K. Schmidt (2006), in prepar 

ration to submit. 
[22] L. Hostlcr and R. Pratt, Phys. Rov. Lctt. 10, 469 (1963). 
[23] L. Hostler, J. Math. Phys. 11, 2966 (1970). 
[24] M. F. Herman, E. J. Bruskin, and B. J. Berne, J. Chem. 

Phys. 76, 5150 (1982). 
[25] G. D. Mahan, M any-P article Physics (Kluwer Aca- 

domic/Plonum Publishcrs, 2000). 
[26] S. Doniach and E. H. Sondheimer, Green's Functions for 



Solid State Physicists (W. A. Benjamin, Inc., 1974). 
[27] W. Kolos and L. Wokiicwicz, J. Chem. Phys. 43, 2429 
(1965). 

[28] K. P. Huber and G. Herzberg, Molecular Spectra and 
Molecular Structure, vol. 4 Constants of Diatomic 
Molécules (Van Nostrand Reinhold, New York, 1979). 

[29] L. Wolniewicz, J. Chem. Phys. 99, 1851 (1993). 

[30] E. Ishiguro, T. Arai, M. Mizushima, and M. Kotani, 
Proc. Phys. Soc. A 65, 178 (1952). 

[31] W. Kolos and L. Wolniewicz, J. Chem. Phys. 46, 1426 
(1967). 

[32] J. Rychlewski, Mol. Phy. 41, 833 (1980). 
[33] P. A. Hyams, J. Gerratt, D. L. Cooper, and M. Raimondi, 
J. Chem. Phys. 100, 4417 (1994). 



